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In this paper, we use a variety of mathematical techniques to explore existence, local stabil- 
ity, and global stability of equilibria in abstract models of mitochondrial metabolism. The 
class of models constructed is defined by the biological description of the system, with min- 
imal mathematical assumptions. The key features are an electron transport chain coupled to 
r^p . a process of charge translocation across a membrane. In the absence of charge translocation 

these models have previously been shown to behave in a very simple manner with a single, 
globally stable equilibrium. We show that with charge translocation the conclusion about a 
unique equilibrium remains true, but local and global stability do not necessarily follow. In 
sufficiently low dimensions - i.e. for short electron transport chains - it is possible to make 
qq ; claims about local and global stability of the equilibrium. On the other hand, for longer 

\0> • chains, these general claims are no longer valid. Some particular conditions which ensure 

stability of the equilibrium for chains of arbitrary length are presented. 

o . 
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1 Introduction 



The processes of electron transport and oxidative phosphorylation in mitochondria 
are of vital biological importance, being central to cellular respiration and hence 
energy production in most eukaryotic cells. Summaries of these processes can be 
found in many modern biochemistry textbooks such as [1] or [2]. The basic features 
of mitochondrial electron transport and oxidative phosphorylation are now well 
understood, but elucidation of many of the detailed mechanisms is still in progress 
[3]. 

1 Funded by an EPSRC/MRC grant to the MIAS IRC (Grant Ref: GR/N14248/01) 
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Mitochondrial electron transport occurs via a series of coupled redox reactions in 
the mitochondrial inner membrane. After the initial reduction of a first electron 
donor (e.g. NADH or FADH 2 produced by glycolysis and the TCA cycle) elec- 
trons are transferred from substrate to substrate, finally being accepted by oxygen. 
During some of these electron transfers a second process takes place - protons 
are pumped across the mitochondrial inner membrane producing a proton gradient 
across this membrane. These protons then return down their gradient, either pas- 
sively (termed leak current) or through a particular enzyme, ATP synthase, leading 
to the phosphorylation of ADR 

Generic models of electron transport chains were explored in [4], where the main 
emphasis was on the input-output response of such models. In the simplest case, 
where the proton gradient across the membrane was ignored, these models were 
found to have very simple behaviour - at all physically meaningful parameter val- 
ues there was a single, globally stable, equilibrium. In [5], this result was shown 
to generalise to the case of electron transfer networks with more general topology 
than a chain. On the other hand in the more biologically realistic case - where the 
build up of a proton gradient has an inhibitory effect on electron transport - analysis 
of the models proved harder. In this paper we analyse in more detail the behaviour 
in this case. 

Before discussing generic models, it is worth mentioning that there are several de- 
tailed models of electron transport and oxidative phosphorylation such as [6], [7], 
[8], [9]. These ordinary differential equation models have been designed with nu- 
merical data in mind, and reflecting the complexity of the processes involved, the 
functional forms are quite involved. Our interest in mitochondria was originally 
inspired by analysis and simulation of some of these numerical models, but the ap- 
proach here is quite different, and more akin to work in [4], [5], [10]. The generic 
model we construct could be instantiated in a great variety of numerical models, 
and the claims we make are valid for all possible instances of the generic model. 



2 The model 



2.1 The basic reaction scheme 



The basic reaction scheme of interest here was described in some detail in [4] but 
will be summarised here. Assume that there are n substrates, each of which can 
exist in an oxidised state Aj and a reduced state Bj so that 

A; + e~ <=> Bi 
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Further, assume that protons can exist in two compartments - the mitochondrial 
matrix (where they are termed H+), and the intermembrane space (where they are 
termed H+) - with the possibility of transfers of the form 

We are interested in reactions which are in general the combination of three pro- 
cesses, a reduction, an oxidation, and the transport of some protons across the mem- 
brane. So for example, if substrate Aj is reduced to Bj, Bj is oxidised to Aj, and p 
protons are pumped across the mitochondrial membrane we get the half reactions 

Ai + e~ ^ Bj, Bj ^ Aj + e~ and pYf m ^ pH + e 
which combine to give 

A + Bj+pH^ A j+ B i + j pH e + 

We also allow the possibility that a reducing/oxidising agent may be external to the 
model giving reactions such as 

A i + pH;^B i + pFC or B ; + pH+ ^ A, + pH + e 

A set of reactions of the kind just described can be combined into a network of re- 
actions. A chain structure (as opposed to a more general network) derives from the 
assumption that each oxidised substrate accepts an electron from only one donor, 
and each reduced substrate transfers its electron to only one acceptor. This intro- 
duces a natural ordering on the substrates, so that for i < n, the z'th substrate is able 
to donate electrons to the (z + l)th substrate, while for i > 1, the z'th substrate is 
able to accept electrons from the (z - l)th substrate. The first substrate is able to 
accept electrons from outside the chain (reflecting the initial reduction of NADH 
or FADH 2 ), and the nth substrate is able to donate electrons to an acceptor outside 
the chain (reflecting the action of 2 ). 

Thus there are n + 1 redox reactions and the z'th reaction has forward rate fi. We 
make no assumptions about the sign of the fi, potentially allowing reactions to be 
reversible. For i < n, the z'th reaction involves the reduction A i5 and for i > 2, the 
z'th reaction involves the oxidation of Bi_! . We define p { as the number of protons 
pumped across the mitochondrial membrane by the z'th reaction. Assuming that the 
quantities p t are constant discounts the possibility of "redox slip" [11], which does 
not appear to be very important in normal circumstances [12]. A quantity \p can 
be defined so that transfer of a single proton across the membrane creates one unit 
of if/, if/ can take any real value and is a strictly increasing function of the elec- 
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trical/chemical gradient against which protons are pumped across the membrane, 
generally termed the proton motive force. 



Finally, reflecting the combined effect of proton leak and ADP phosphorylation, 
there is a process with rate L representing the "decay" of if/. When there is no 
gradient, no protons leak through the membrane, so that L(0) = 0. Further L is 
assumed to be strictly increasing in if/. 

The structure of the model is illustrated in Figure 1 . 




n+l 



Fig. 1. A schematic representation of the reaction network. The quantities A; and Bj refer 
to oxidised and reduced states of the substrates. The functions /; define the forward rates 
of reaction of the n + 1 coupled redox reactions. The quantity if/ represents the electrical 
and chemical gradient across the mitochondrial membrane, which has an inhibitory effect 
on any redox reactions which involve proton pumping. 

Because the total quantity - oxidised plus reduced - of any substrate in the chain 
is conserved, reduced forms of the substrates are not explicitly introduced. Instead, 
the concentration of Aj is referred to as x t , and the total concentration of A + Bi is 
assumed constant at m,. We arrive at a model of the form: 



X\ = -fl(x U l/f) + f2(Xi,X 2 ,lff) 

Xi = -fi(Xi-i, x h iff) + f i+l (x h x i+u iff) i = 2,...,n-\ 

Xn fn(x n —\,X n , 
n+l 

* = 2 Pifi - U&) 

i=l 



(1) 



The phase space of this system is defined by the equations: 



< Xi < nit i = 1 , . . . , n 
— oo < if/ < oo 

and is hence n + l dimensional, being the product of a closed n-dimensional box 
and the real line. 
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2.2 Assumptions 

All the functions f h along with L, are assumed to be C 1 (once differentiable in all 
their arguments with continuous derivatives). The following notation is used for the 
derivatives of the functions fi: 

_ dfi _ _ dfi 

fij = ' Fij = ~fij ' fill' = ' = ~f<i(> @) 

At finite substrate concentrations, all reaction rates are finite, so that at any fixed if/ 
each is bounded on its domain of definition. 

Since ifj represents a potential against which some of the reactions must do work, 
the following relations are obtained: 

^<0ifp«*0 and /^ = 0if/>j = (3) 

If pi ± 0, then if/ inhibits the forward reaction and we assume that sufficiently large 
values of if/ make the reaction rate arbitrarily small or negative, i.e. 

lim fi(;if/) < i = l,n + 1 
lim fi(-,; if/) < i = 2,...,n 

This reflects the fact that the energy required to pump a proton against a chemi- 
cal and electrical gradient becomes large as the gradient increases. Similarly -if/ 
inhibits the backward reaction so that: 

lim fi-, <A) >0 i = l,n + l 

ft— >— CO 

lim fi(;;ffr) > i = 2,...,n 

The following equations imply that no reaction can proceed in the absence of any 
of its substrates: 

/i(0,-) = 
M;0,-) = i 
fiinti-u-,-) = i 
f n+ \(m n ,-) = 



2, •••,« 
2, •••,« 



(4) 
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The final set of conditions imply that increased substrate concentration increases 
the rate of reaction unless one of the substrates is entirely absent: 



fu >o 

f u > and f u > if x M < m^ x i = 2, • • • , n 
fi+u < and f +u > if x i+l > i = 1, 1 



(5) 



The fact that the first and final inequalities are always strict implies that there is al- 
ways some electron donor to reduce the initial substrate, and some electron acceptor 
to oxidise the final substrate, and ensures nondegenerate behaviour. The assump- 
tions from (5) mean that f u , Fjj and as defined in (2) are all nonnegative. The 
definition of these nonnegative quantities is solely to simplify later arguments. 



3 General behaviour of the system 



In this section we outline some properties of the model that hold regardless of the 
number n of redox pairs. 



3.1 Boundedness of solutions 



It is convenient to define an n x (n + 1) matrix which can be regarded as a stoichio- 
metric matrix for the redox reactions: 



S = 



-1 1 
-1 








1 1 



Defining the vector of reactant concentrations x = [xi , x 2 , ■ ■ ■ , x n ] T , the vector of re- 
action rates v(x, if/) = [fufo, . . .f n +i] T , and the nonnegative vector P = [p u . . . ,p n +i~\ T , 
we can rewrite the system of equations (1) more briefly as 



x = S v(x, if/) 

(A = J P r v(x,<A)-L(iA) 
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We now show that all forward trajectories of the system are bounded. Since the 
phase space is bounded in x, what needs to be shown is that all trajectories en- 
ter a bounded region in the if/ direction. This amounts to showing that ifj > 
for if/ sufficiently large and negative, and that if/ < for if/ sufficiently large and 
positive. By assumption, for any given i, either p t = or is strictly nega- 
tive and lim^oo /;■(•, -,^) < 0, lim^-a, _/,(-, -,if/) > 0. This in turn implies that 
lim^_>oo P T \(x, if/) < and lim^o, P T \(x, if/) > 0. In addition L is strictly increasing 
from zero as if/ increases. Thus for any fixed value of x, lirn^oo P T \(x, if/) - L(if/) < 
and lim^-oo P T \(x, if/) - L(if/) > 0. Define if/ (x) as the value of if/ at which 
P r v(x, if/) - L(if/) = 0. iAo(x) is uniquely defined since P T \(x, if/) - L(if/) is strictly 
decreasing. By the implicit function theorem, if/o(x) is a differentiable function since 
P T \(x, if/) - L(if/) is a differentiable function of x. Since it has a compact domain, 
i^o(x) achieves a maximum value which we call ifj max , and a minimum value which 
we call if/ min . By these definitions, if/(if/, x) < for all if/ > if/ max , and ij/(if/, x) > for 
all ifj < ifj min . 

Thus all trajectories enter a closed box, S, bounded by the hyperplanes = 0, 
xi = mu if/ = ifj min and if/ = if/max, and this box forms a trapping region for the system 
in all dimensions. 



3.2 The Jacobian 



Direct calculation gives that the Jacobian, J, of the system is: 



J = 



-fxx ~ Fix fn 

Fix -fn - F32 



Pxfxx-PiFix Pifn-P3F 32 








Fxi)/ ~ Fjjf, 
Fiijj ~ F34, 



fnn F n +\ n 



Fnip F n+ \fp 
n+l 

Pnfnn ~ Pn+ X F n+ \ >n — ^ PiFty 

i=X 



Here = ^. The structure of this Jacobian can be made clearer by defining 
two further quantities: A nonnegative vector in M", F = [F^, . . .,F ni// ] T ; and an 
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(n + \)xn matrix 



v ~ J i- 

ox 



fn 





•• 


• 


-F 21 




•• 


• 





-F32 


/33 ' ' 


• 








•• 


fnn 








•• 


~Fn+\,n 



Then the Jacobian can be written in the block form: 



7 = 



SV SF 

pTy _pTf _ T 



(6) 



S V is the Jacobian of the system without feedback, which is tridiagonal, and can 
easily be shown to have real negative eigenvalues [4]. It was shown in [13] that the 
structures of S and V along with the nonnegativity of P and F imply that 7 is a so 
called / J(_) matrix (see Appendix A for the definition!!]. This result is independent 
of n, the length of the chain. It has the consequence that the system is injective; this 
is discussed further in the next section. 

The fact that 7 is a i ,< ~ ) matrix has another consequence of importance to us: It 
means that its eigenvalues are excluded from a certain wedge around the positive 
real axis: If A = re' e is an eigenvalue of an m x m P matrix, then it is proved in [14] 
that: 

\6 — tc\ > n/m 



and equivalently for a 7 ,( } matrix, 



\6\ > n/m 



Clearly when m = 2, this means that both eigenvalues lie in the left half plane, 
so that 2x2 7 J( " ) matrices are Hurwitz stable (see Appendix A for a definition of 
"Hurwitz stable" which we will abbreviate to "Hurwitz"). However for m > 2, 7 J(_) 
matrices may be unstable. 



2 The nondegeneracy conditions presented in [13] are met because the nth substrate is 
terminal, and all substrates are able to transfer electrons along the chain to the nth substrate. 



8 



3.3 A unique equilibrium 

The existence of a unique equilibrium for this system was shown in [4] by a direct 
method. It also follows from the arguments presented above: That an equilibrium 
must exist follows, by the Brouwer fixed point theorem, from the existence of the 
compact, convex, trapping region, S constructed above; That this equilibrium must 
be unique follows from the fact that the Jacobian is a matrix, and hence the 
system is injective [15]. Thus as our first result we can state that 

Electron transport chains coupled to charge translocation across a membrane 
have exactly one equilibrium. 

It is interesting that the possibility of multi stability is immediately ruled out. How- 
ever this in itself does not tell us whether all trajectories must necessarily converge 
to the unique equilibrium, or whether periodic or chaotic behaviour is still possible. 



4 Stability of the equilibrium 

In this section, we analyse stability of the unique equilibrium, starting with low 
dimensions (i.e. short chains). For two dimensions we prove that the equilibrium 
is globally asymptotically stable. In three dimensions we show that the addition 
of an extra, reasonable, constraint implies that the equilibrium is locally stable, 
and further constraints ensure that it is globally stable. We then demonstrate that 
these constraints do not suffice to guarantee stability in four dimensions and higher. 
Finally, we outline some additional special conditions that guarantee the Jacobian 
is Hurwitz in all dimensions. 



4. 1 The system in two dimensions 

The system in 2D consists of a single redox pair subject to a reduction process 
and an oxidation process, both possibly coupled to proton translocation across the 
membrane. It takes the form 

x 1 = -f 1 (x 1 ,if/)+f 2 (x 1 ,iff) 
<A = Pi/i + Pifi - U^) 

The Jacobian of the system in this case is: 
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h = 



-fn ~ Fix - F 2l /, 

P\f\\-piF 2 \ -L^-p\F\,p -p 2 F 2l p 



(7) 



We have already mentioned that 2D matrices are Hurwitz stable, and it fol- 
lows that the matrices J 2 are Hurwitz stable (This can also be shown with a direct 
calculation). 

Since J 2 is Hurwitz stable everywhere, not just at the unique equilibrium, the Markus- 
Yamabe Theorem (e.g. [16], [17], [18]) ensures that the equilibrium is globally 
stable. We also offer an alternative, elementary, proof of global stability. By the 
Poincare-Bendixson Theorem (see, for example, [19]), w-limit sets of a flow on 
compact subsets of R 2 must either contain equilibria or consist of a periodic orbit. 
In this case we can rule out the possibility of periodic orbits: The divergence of the 
vector field is equal to 



which is negative. Thus the vector field satisfies the Dulac criterion (e.g. [20]) and 
there are no periodic orbits. We know that there is only one equilibrium, which is 
locally stable, and therefore there are no heteroclinic or homoclinic orbits either. 
Since every forward trajectory enters the box S, the unique equilibrium must be 
the to- limit of every trajectory, and is hence globally stable. 

4.2 The system in three dimensions 

Slightly more complex than the two dimensional system is the system in three 
dimensions which takes the form 

x\ = -fi (xi , if/) + f 2 (Xi ,x 2 ,ifj) 

X2 = -fliXi, X 2 , (A) + MX2, <A) 

<A = Pif\ + Pifi + Psfs - L(if/) 
with Jacobian 



Tr(J) = -fn- F 2l -p\F^- p 2 F 2ip - 



~fu - F 2 \ f 22 Fty - F 2l// 

F 21 —fn - F 32 F 2l// - F 3lf/ 

P\f\\-piF 2 \ Pifn-p-iF^ -L^-p\F\^ -p 2 F 2lf/ -p^F^ 



(8) 
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As it stands, 73 is not always Hurwitz. For example, the Jacobian constructed using 
the following values: p x =3,p 2 = 0,p 3 = 88, F llf , = 33, F 2l/ , = 4,F 3lf , = 0.6, f n = 
23, f 2 2 = 3, F 2 \ = 94, F 32 = 76, = 6 has two eigenvalues with positive real part. 



7 3 can be shown to be Hurwitz everywhere in 3D provided one extra condition is 
met: p x and p 3 must have the same ordering as F^ and F 3l// . For a real number z, 
define the function 



sign(z) = 



1 (z>0) 
(z = 0) 
-1 (z<0) 



(9) 



Then the ordering assumption translates to the following statement: 

sign(7V - F ii) = sign(p 3 - Pi) (10) 



With this assumption, the Jacobian is everywhere Hurwitz, and hence the equilib- 
rium is locally asymptotically stable. The proof is simple but requires some lengthy 
evaluations, and the details are presented in Appendix B. 

Unlike in the 2D case it does not follow that the equilibrium is globally stable, since 
the Markus-Yamabe conjecture does not hold in dimensions greater than 2 [21]. 
However we can prove global stability in this case too subject to a strengthened 
version of the ordering assumption on the quantities p t and F^. We now require 

sign(/> - F # ) = signO -pj) (11) 



far 1, ye {1,2,3}. 

With this assumption we are able to use a version of Li and Muldowney's au- 
tonomous convergence theorem (Theorem 4.1 in [22]) to show that the unique 
equilibrium is globally stable. In order to use this theorem two concepts are needed: 

(1) The second additive compound of a matrix 

(2) Logarithmic norms of a matrix 

Both quantities are defined for square matrices. The second additive compound ma- 
trix of any nx n matrix 7 is a square matrix of dimension n C 2 which we will term 
7 [2] . Logarithmic norms are scalar quantities, and corresponding to any given ma- 
trix norm, there is a logarithmic norm. Unlike matrix norms, however, logarithmic 
norms may take negative values. The definitions are given in Appendix A. 
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Consider a dynamical system with Jacobian J(x) at some point of phase space x. 
Define J to be the set of all these Jacobians. For our purposes, the autonomous 
convergence theorem states the following: If a logarithmic norm ji can be found 
such that 



H(J [2] ) < for all J e J 



(12) 



then the limit set of each bounded semi-trajectory of the dynamical system is an 
equilibrium. 

Since all trajectories enter the trapping region & in our system, and since S contains 
a unique equilibrium, finding a suitable logarithmic norm satisfying (12) will suffice 
to prove global stability of the equilibrium. 

The second additive compound in this case is: 



r[2] _ 
J 3 ~ 



-f\\-F2\ -fn-F, 



32 



Pihi-p^F^i -fu-Fn -F^-Ya PiFty 



i+i 



-(Pifn-PiFii) 



-(Fty - F^) 

hi 



'21 



-fn ~ F32 -Lf - H PiFty 
(+1 



We will construct a logarithmic norm /u T such that [i T (jf^) < 0. For a real n x n 
matrix, the logarithmic norm corresponding the usual || • ||i norm takes the form: 



Hi = max 

ie{\,...,n] 



Xu + 



2 \ x ki 



From the definition it is clear that a matrix has negative logarithmic norm //1 if and 
only if every diagonal entry is negative and it is strictly diagonally dominant in 
every column. Next we define a constant diagonal coordinate transformation 



T = 



10 



0^0 

Pmax 

0^ 

Pmax 



where p max = max (pi). 

(€{1,2,3} 

According to Lemma 2.2 of [23], given any invertible transformation T, jJ. T (M) = 
Pi(TMT~ { ) defines a new logarithmic norm. In this case, since T is a diagonal 
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prove that /u T (J 1 ^) < 0, we need to show that J' = TJ [ ^ ] T 1 is strictly diagonally 



matrix, the diagonal entries of M are the same as those of TMT 1 . Thus in order to 
prove that /u T (jf ] ) < 0, we 
dominant in every column. 

For the first column, we have 



Ai + K21I + K31I - ~fi2 ~ F32 - fn - F 2 i 

+ 



Pi , P3 
J22 ^32 



+ 



P2 Pi , 
^21 ill 



It can easily be seen that the term on the right hand side is negative since for any 
two nonnegative scalars \a-b\< max{|a|, \b\\. 

For the second column, we have 

3 

J' 22 + \ j' l2 \ + \ j' 32 \ = - ^^PiFty -1,^- f n + p max \F 2i /, - F 3 ^\ 



i=i 



For the final column, we have 



^33 + |^13| + 1^23 1 - ~^jPiFi^ - - F 32 + Pmax ~ ^li/r| 



In order to show that the right hand sides of the last two expressions are negative 
we need to show in each case that our ordering assumption (11) implies that the 
final term (which may be positive) is dominated in magnitude by the other terms. 

Note that W^-F ;A < maxjF,-^, F^} < max (F w ). Then there are only three cases: 

/te{l,2,3| 



(1) if Pmax 


= P\, then p max 


Fit/, 


- Fty 


< p\F H , and p max 


Fit/, 


- F^ 


< P\F\^ 


(2) if/w 


= p 2 , then p max 




~ F 3l p 


< PiF^, and p max 


Fity 


- Fi^ 


< PiFty 


(3) ifp max 


= p 3 , then p max 


Fitp 


- F^ 


< p 3 F 3l// , and p max 


Fiif, 


- Fty 


< PiF^ 



Each of these possibilities leads to the same conclusion - that J' u + £ l^jyl < f° r 



each i. Hence we have fi T yjf ] j < 0. 



k,k+i 



This result means that if the ordering assumption (11) holds, then the unique equi- 
librium is globally stable. The ordering assumption itself has the following reason- 
able physical meaning which we would expect to be fulfilled in practice: If redox 
reaction i is involved in pumping more protons across the membrane than redox 
reaction j, then reaction i is correspondingly more inhibited by \p than reaction j. It 
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is interesting to note however that this assumption is not necessary to prove global 
stability in the 2D case. It is also unknown to us whether the weaker assumption 
(10), which guarantees that the Jacobian is everywhere Hurwitz, actually guaran- 
tees global stability in 3D. 

4.3 Unstable examples in higher dimensions 

The ordering assumption (11) does not guarantee global or even local stability of 
the equilibrium in dimensions greater than 3. It is easy to construct counterex- 
amples. For example, in four dimensions, the Jacobian constructed by choosing 
P\ = 2, p 2 = p 3 = 0, p 4 = 73, = 167, F 2l p = F^ = 0, F^ = 176, f n = 4, 
fn = 7, / 33 = 1, F 2 i = 32, F32 = 64, F 43 = 174, = 33, satisfies all the con- 
straints, including the ordering assumption on the values of p t and F ilf ,. However it 
has, two eigenvalues with positive real part. 

We make the following remarks: 

(1) By continuity, the fact that a non-Hurwitz Jacobian can be constructed in 4 
dimensions guarantees that such examples also exist in all higher dimensions. 

(2) Systems with non-Hurwitz Jacobian satisfying the ordering assumption (11) 
seem to be rare. Through use of an automated computer script running in the 
open source numerical computation program Scilab [24], counterexamples in 
dimension 4 were found by randomly choosing values for the different terms 
in the Jacobian, such that all the assumptions were satisfied. Out of hundreds 
of millions of sets of values, less than ten were non-Hurwitz. 

(3) The counterexamples found appear always to be close to breaking the ordering 
assumption. For instance, in the example shown, p A is much greater than p x , 
whereas F^ is close in magnitude to F^. 

4.4 A special case: Reaction rates dependent on potentials 

In this section we consider an interesting assumption which ensures that the Jaco- 
bian is Hurwitz everywhere (and hence the unique equilibrium is locally stable). 
The assumption is as follows: 

(1) Associated with each half reaction is some "potential": In the case of a redox 
reaction of the form Aj + e~ ±=? B i? a potential means any strictly increasing 
scalar function of [A;]; In the case of a charge transfer across a membrane a 
potential means any strictly increasing scalar function of ijj. 

(2) The rate of any full reaction depends only on the sum of the potentials for the 
half reactions involved, and is a strictly decreasing function of this sum. 
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This assumption can be interpreted, loosely, as saying that the energetics of the 
system determine the reaction rates. For example, consider the electron transfer 
coupled to some proton pumping 



Ai + Bj+pH^Aj+Bi+pH: 



derived from the half reactions 



A + e Bi, Bj ^ Aj + e and pYt m <=; pR + e 



In this case, the assumption would imply that the forward rate of the combined 
reaction can be written f(-gj(xj) + gi(Xi) - pg^iif/)) where the only stipulation is 
that /, gi, gj and g^ are strictly increasing in their arguments. When this assumption 
is made about all reaction rates in the system, the full system becomes: 



xi = -/i(giOi) - Pig^)) + fi{-g\(x\) + g 2 (x 2 ) - Pig^)) 

Xi = -fi(~gi-l(Xi-l) + gi(Xi) - Pig^)) + 

fi+i(-gi(xd + gi+ifei) - pMg^)) i = 2,...,n 

X„ = -fn(-gn-l(Xn-l) + gn(x„) ~ Png^)) + fn+\(-gn(x n ) ~ Pn+\g^)) 
n+l 

^ = J]p i f i -L(^ 

(=1 



The term /v(-g,-i(X-i) + gi(xd - Pig^)) represents the rate at which the z'th sub- 
strate receives electrons from the (z - l)th substrate. Denoting by f., g\ and g, the 
derivatives of the functions f h g t and g^, the Jacobian of this system can be written 
J = JqD where J is the symmetric matrix 



fi -(/ 2 +/ 3 ') 



Pxf'x-Pifi 








P\f[- Pif 2 
Pif' 2 ~ Pifs 



-(fn + f n +l> Pnfn ~ Pn+lf, 



n+l 



Pnfn~ Pn+lf „ +1 ~Z P]fi ~ ~f 
■ 1 



(13) 
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and D is the positive diagonal matrix 




••0 



••0 



D = 



(14) 











■■g'n 



From the discussions earlier, J$ is a matrix. Further it is symmetric, and hence 
sign symmetric (see Appendix A for a definition of sign symmetry). This implies 
[25] that J is D-stable, i.e. the product of J with any positive diagonal matrix is 
Hurwitz. Hence J is Hurwitz. Thus the assumption that reaction rates depend on 
the sum of potentials of the half reactions involved ensures that the Jacobian of the 
system is everywhere Hurwitz. 



5 Discussion and conclusions 



We have analysed in some detail, and using a variety of mathematical techniques, 
the behaviour of electron transport chains coupled to a charge translocation process. 
In all cases trajectories are bounded, and there is a unique equilibrium, but ques- 
tions about the stability of this equilibrium have proved harder. Where the chain 
consists of a single redox pair, the unique equilibrium is globally stable. When 
there are two redox pairs the same conclusions can be reached subject to some 
extra conditions on the feedback process. In higher dimensions no such general 
conditions could easily be found. Thus the length of the electron transport chain is 
crucial in deciding on stability of the equilibrium. 

It is somewhat surprising that the coupling of electron transfer to a membrane po- 
tential - a negative feedback loop - can serve to destabilise the unique equilibrium 
in these systems. Interestingly, when the reaction rates are monotonic functions of 
a sum of potentials, then the system in any dimension could be proved to be ev- 
erywhere Hurwitz. Reaction rates cannot in general be seen in this way, but in the 
case of reactions which are primarily about charge transfer, the assumption could 
be reasonable. Certainly some of the choices of reaction rates in numerical models 
such as [6] satisfy this assumption. 

There are some interesting open questions, both biological and mathematical. From 
a biological point of view, it is of interest to find out whether experiments on mi- 
tochondria with constant inputs ever display behaviour other than convergence to 
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an equilibrium, such as periodic or chaotic behaviour. If this is never the case, then 
this suggests that our very general model may be omitting certain important biolog- 
ical/thermodynamic restrictions on the reaction rates, which would tend to stabilise 
the system. It would also be interesting to see how additional processes such as 
transport processes in the full numerical models ([6], [9] for example) affect the 
conclusions presented here. 

An open mathematical question is whether there are equivalent conditions to the 
ordering condition in 3D which ensure that the Jacobian of the system is Hurwitz 
in arbitrary dimension, or better still that the second additive compound has nega- 
tive logarithmic norm, and hence the unique equilibrium is globally stable. If such 
conditions exist can they be given general biological meanings? 

It would also be interesting to explore when the results presented here survive 
weakening of the assumption that electrons are transferred along a chain. Although 
electron transfers taking place in the mitochondrial membrane are often described 
via a "chain" it is likely that this description is a convenient simplification rather 
than the whole truth. General electron transfer networks in the absence of a poten- 
tial were analysed in [5] and found to have simple behaviour. Application of the 
theory presented in [13] should allow determination of when these networks give 
rise to Jacobians when interacting with a membrane potential. 

Finally, although conditions ensuring sign-symmetry of the system imply that the 
Jacobian is everywhere Hurwitz, it is an open question as to whether this implies 
global stability of the unique equilibrium. Since the Markus-Yamabe conjecture 
does not hold in dimensions greater than 2 [21], global stability does not follow 
automatically from local stability, and requires independent proof. 



A Definitions 



A. 1 Hurwitz stability of matrices 

A square matrix is defined to be Hurwitz stable if all its eigenvalues lie in the open 
left half of the complex plane - i.e. the real parts of all its eigenvalues are negative. 



A.2 P matrices and related classes 

For some nxm matrix A, A{a\y) will refer to the submatrix of A with rows indexed 
by the set a c { 1 , . . . , n\ and columns indexed by the set y c { 1 , . . . , m}. A principal 
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submatrix of A is a submatrix containing columns and rows from the same index 
set, i.e. of the form A(a\a). A minor is the determinant of any square submatrix 
of A. If A{a\y) is a square submatrix of A (i.e. \a\ = |y|), then A[ar|y] will refer to 
the corresponding minor, i.e. A[ar|y] = det(A(a|y)). A principal minor of A is the 
determinant of a principal submatrix of A. 

P matrices are square matrices all of whose principal minors are positive. They are 
by definition nonsingular. If -A is a P matrix, then we will say that A is a 
matrix. If A is a matrix, this means that each k x k principal minor of A has sign 
(-!)"• 



A3 Sign symmetry 



An n x n matrix is sign-symmetric if symmetrically placed minors have the same 
sign, i.e. A[a|y]A[y|Q'] > for every a,y C {1, . . . ,n] with \a\ = |y|. 



A.4 Second additive compound matrices 



A brief definition of the second additive compound of any square matrix can be 
found in [26]. For a more detailed discussion see [27]. For a 3D matrix 



( \ 
a\\ a\i an 

a%\ a22 o-n 

«31 £32 #33 



(A.l) 



the second additive compound takes the forn{EI 



l PI 



an + ^22 «23 — ^13 

^32 ^11 "I" ^33 ^12 
—a^\ an «22 + ^33 



This second additive compound was constructed using the standard lexicographic 
ordering of basis vectors. It is possible to construct a second additive compound 
using a different ordering, but such choices make no difference to the logarithmic 
norms of the matrix. 

3 In general, the second additive compound of a matrix A has dimension d C2 where d = 
dim(A). When dim(A) = 3, we get dim(A [2] ) = 3 also, but this is not generally the case. 
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A.5 Logarithmic norms 



If || • || denotes a vector norm on W, and also the induced matrix norm on«xn 
matrices, then the logarithmic norm [28], also known as a Lozinskii measure, of an 
nxn matrix A is defined by 

(A , r \\I + hA\\-l 
H{A) = hm (A.2) 

fc->0 + h 



B Local stability in 3D 



In this appendix we prove local stability of the equilibrium in three dimensions, 
subject to the assumption in (10), using the Routh-Hurwitz theorem. Consider the 
characteristic polynomial of a matrix A: 

\AI-A\ = A" + M"" 1 + ■ ■ ■ + b n -iA + b n (B.l) 



In this equation, / is the n x n identity matrix, and the coefficients Z? ( are the sums 
of all principal minors of -A of dimension /. For a matrix, bt > for all i. Now 
define bk = for all k > n, and construct a set of numbers A, as follows: 



bi 


1 











• 


•• 


h 


b 2 


bi 


1 





• 


•• 


b 5 


b 4 


b 3 


b 2 


by 


1 • 


•• 



: : : "-. : 

bn-\ bn-i b 2 i-3 b 2 i-4 b 2i - 5 b 2i -e ■ ■ ■ b t 



The Routh-Hurwitz theorem states that A is Hurwitz if and only if A, > for all 
i < n. In three dimensions, we need to check that the three quantities 



A l= b, (B.3) 
A 2 = b x b 2 - b 3 (B.4) 
A 3 = b 3 (b l b 2 -b 3 ) = b 3 A 2 (B.5) 

are all positive. Since all the b t are positive, all three quantities are positive if and 
only if A 2 > 0. This in turn follows (condition 12 in [25]) if 

< ai 2 a 23 a 3 i + a 2 \a 32 a\ 3 - 2a\\a 22 a 33 
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where a, 7 are elements of A. Substituting a, 7 for the elements of the Jacobian and 
expanding using the open source symbolic algebra program Maxima [29] gives the 
following condition: 



«12«23«31 + «21«32«13 _ 2ana 2 2«33 - F 2 l ^32 \2p 3 F 3l p + ZpiFty ~ P^Xxf,) 

/22 (2^3^ + ^P\F^ - P\F 3ip ) 



+ positive terms 



With the ordering assumption (10), we get: 



2p 3 F 3 ^ + 2p x F^ - p 3 F lip > 
2p 3 F 3lfr + 2p\F Xil - p x F 3ljj > 



(B.6) 
(B.7) 



Thus the Jacobian is everywhere Hurwitz and hence the unique equilibrium of 
the system must be locally asymptotically stable. Note that the restriction (10) is 
stronger than necessary to ensure that J is Hurwitz, but no other set of conditions 
with a clear physical meaning that make the Jacobian Hurwitz have been discov- 
ered. Finding a set of necessary and sufficient conditions for J to be Hurwitz is a 
difficult problem. 
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